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1.  INTRODUCTION 


The  field  of  detonation  physics  has  enjoyed  considerable  attention  since  the  phenomenon  was 
recognized  over  a  century  ago  (Berthelot  and  Vielle  1881, 1882a,  1882b;  Mallard  and  LeChatelier  1881). 
Benefits  that  would  result  from  the  ability  to  exploit  and  control  such  energetic  events  have  led  to 
numerous  explorations  designed  to  determine  the  character  of  a  system  and  conditions  that  will  result  in 
a  detonation.  The  extreme  pressures  and  high  energies  released  in  a  detonation  have  posed  considerable 
experimental  challenges  toward  these  characterizations.  Even  greater  obstacles  in  monitoring  detailed 
microscopic  chemical  and  physical  changes  in  the  detonation  are  the  time  and  length  scales  over  which 
the  event  occurs.  These  experimental  challenges  have  necessitated  development  of  theories  that  attempt 
to  describe  the  phenomenon  and  complement  the  experimental  analyses.  The  majority  of  theoretical 
treatments  are  based  on  hydrodynamic  theories  that  predict  the  macroscopic  behavior  of  a  system  that 
detonates,  such  as  changes  in  pressure,  temperature,  and  density  (Fickett  and  Davis  1979;  Fickett  1985). 
These  models,  however,  do  not  give  information  about  the  microscopic  chemical  and  physical  processes 
occurring  during  this  violent  event  and  in  many  instances,  rely  on  substantial  approximation.  Very  little 
is  known  about  the  chemical  reactions  that  initiate  and  sustain  a  detonation,  although  theoretical  (Tsai  and 
Trevino  1984;  Tsai  1990;  Blais  and  Stine  1990;  Karo,  Hardy,  and  Walker  1978;  Elert  et  al.  1989; 
Kawakatsu  and  Ueda  1988, 1989;  Kawakatsu,  Matsuda,  and  Ueda  1988;  Brenner  1992;  Robertson  et  al. 
1992;  White  et  al.  1994;  Brenner  et  al.  1993)  and  experimental  methods  (Gupta  et  al.  1995;  Fried  and 
Ruggiero  1994;  Dlott  and  Fayer  1990;  Tokmakoff,  Fayer,  and  Dlott  1993;  Chen,  Tolbert,  and  Dlott  1994; 
Chen  et  al.  1995;  Hong,  Chen,  and  Dlott  1995)  are  being  directed  toward  unraveling  details  of  these 
ultrafast,  violent  events. 

The  method  of  molecular  dynamics  is  a  powerful  and  well-established  simulation  technique  to  provide 
details  of  the  atomistic  processes  occurring  in  a  chemical  reaction  (Raff  and  Thompson  1985).  Though 
the  majority  of  molecular  dynamics  simulations  have  focused  on  gas-phase  problems,  computer  power  has 
increased  to  the  point  that  molecular  dynamics  simulations  of  condensed  phase  can  be  realized.  As  early 
as  the  1970s,  molecular  dynamics  simulations  of  simple  solid  models  predicted  reasonable  features 
expected  from  a  shocked  solid  (Holian  et  al.  1980;  Powell  and  Batteh  1979,  1980).  More  sophisticated 
models  incorporated  energy  release  reactions  that  could  predict  detonation  (Tsai  and  Trevino  1984;  Tsai 
1990;  Brenner  1992;  Robertson  et  al.  1992;  White  et  al.  1994;  Brenner  et  al.  1993);  however,  the  models 
either  described  the  chemistry  that  drive  the  detonation  qualitatively  incorrectly  (Tsai  and  Trevino  1984; 
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Tsai  1990;  Karo,  Hardy,  and  Walker  1978),  or  the  models  had  undesirable  features  (Brenner  1992; 
Robertson  et  al.  1992;  White  et  al.  1994;  Brenner  et  al.  1993). 

It  is  our  intention  to  use  the  method  of  molecular  dynamics  to  investigate  the  microscopic  processes 
that  occur  during  a  detonation  and  to  determine  the  properties  of  the  system  that  affect  this  phenomenon. 
It  is  also  hoped  that  relevant  mechanisms  will  be  revealed  in  the  process.  Before  we  can  do  this,  we  must 
first  develop  a  model  that  more  correctly  describes  exothermic  chemical  reactions,  and  determine  whether 
it  can  adequately  simulate  the  phenomenon  known  as  a  detonation.  We  have  developed  such  a  model, 
the  features  of  which  are  described  in  the  accompanying  paper  (Rice  et  al.,  in  press).  The  focus  of  the 
study  presented  here  is  to  compare  our  molecular  dynamics  simulations  of  a  shock-initiated  reacting  crystal 
with  predictions  of  classical  hydrodynamic  theory  of  detonation.  Details  and  results  of  the  two  types  of 
calculations  will  be  presented  and  a  comparison  given. 

The  results  of  the  first  calculation,  the  determination  of  the  equation  of  state  of  the  system,  will  be 
used  to  evaluate  the  classical  conservation  equations  that  relate  the  mass,  momentum,  and  energy  of  the 
quiescent  crystal  with  the  state  behind  the  detonation  wave  (Fickett  and  Davis  1979;  Fickett  1985).  The 
three  conservation  equations  are: 

Conservation  of  Mass:  PqD  =  p(D-u)  (1) 

where  p0  is  the  density  of  the  undistuibed  crystal,  p  is  the  density  of  the  system  behind  the  shock  front, 
D  is  the  velocity  of  the  wave  propagating  through  the  undisturbed  crystal  (constant  for  an  unsupported 
detonation),  and  u  is  the  velocity  of  the  products  behind  the  detonation  wave; 

Conservation  of  Momentum:  p  -  p0  =  p0uD,  (2) 

where  p  and  p0  are  the  pressures  of  the  products  behind  the  shock  front  and  quiescent  crystal,  respectively. 
The  variable  u  can  be  eliminated  from  Equations  (1)  and  (2)  resulting  in  the  Rayleigh  line 


R  =  p02D2  -  (p  -  p0)/(v0  -  v)  =  0. 


(3) 
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The  final  equation  requiring  conservation  of  energy  is 


Conservation  of  Energy:  e  +  pv  +  1/2(D  -  u)2  =  e0  +  p0v0  +  1/2  D2  ,  (4) 

where  e,  eo,  and  v,  v0  are  the  specific  internal  energies  and  specific  volumes  of  the  final  and  initial  states, 
respectively.  The  term  "specific"  as  used  here  with  respect  to  some  quantity  refers  to  that  quantity 
normalized  to  unit  mass.  The  variables  D  and  u  can  be  eliminated  from  Equation  (4)  using  Equations  (1) 
and  (2)  to  a  form  that  is  referred  to  as  the  Hugoniot  function. 

h(T,v)  =  e  -  e0  -  l/2(p  -  p0)(v0  -  v).  (5) 

The  set  of  (T,  v)  for  which  Equation  (5)  is  zero  make  up  the  h(TH,  vH)  curve  known  as  the  Hugoniot 
curve.  The  intersection  of  Equations  (3)  and  (5)  determine  the  state  of  a  system  (p,  v)  for  a  given 
detonation  velocity,  D.  There  are  a  series  of  Rayleigh  lines,  defined  by  the  parameter  D,  that  intersect 
the  Hugoniot  curve;  all  but  one  give  two  solutions  to  Equations  (3)  and  (5)  and  represent  unsteady  shocks. 
The  velocity  that  uniquely  satisfies  Equations  (3)  and  (5),  called  the  Chapman-Jouguet  velocity, 
corresponds  to  the  p  -  v  point  where  the  Rayleigh  line  is  tangent  to  the  Hugoniot  curve.  The  C-J  point 
is  the  state  of  the  system  corresponding  to  an  unsupported  detonation,  the  event  we  will  attempt  to 
simulate.  We  will  determine  the  C-J  point  and  detonation  velocity  for  comparison  with  the  molecular 
dynamics  simulation. 

The  second  calculation  is  the  molecular  dynamics  simulation  of  a  plate  impacting  the  quiescent 
diatomic  molecular  crystal,  and  initiating  reaction.  This  will  be  denoted  throughout  as  the  computer 
experiment  We  will  compare  the  predicted  Chapman-Jouguet  point  and  velocity  with  the  shock  wave 
velocity  and  state  of  the  system  from  our  computer  experiment 

This  report  will  first  briefly  describe  the  model,  the  details  of  the  calculations  and  conclude  with 
results,  discussion,  and  comparison  of  the  two  predictions.  All  calculations  described  hereafter  are 
two-dimensional. 
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2.  MODEL 


The  interaction  potential  used  to  describe  the  two-dimensional  crystal  of  diatomic  molecules  arranged 
in  a  herringbone  lattice  is  similar  to  that  used  by  Brenner  et  al.  (1993). 

v  -  E  £  {  W  (2  -  5^)  VR(r„)  -  5T  VA(rJ  ♦  Vvdw  }  (6) 

The  functions  and  parameters  used  in  Equation  (6)  are  given  in  Table  1.  The  features  of  this  potential 
energy  function  are  detailed  in  the  accompanying  paper  (Rice  et  al.,  in  press).  The  low  temperature, 
ambient  pressure  lattice  parameters  for  a  unit  cell  of  this  crystal  in  the  x  and  y  directions  (denoted  a  and  b, 
respectively)  are  4.34  and  6.27  A.  A  unit  cell  is  illustrated  in  the  inset  of  Figure  1.  The  unit  cell  contains 
two  molecules;  center-of-bond  (COB)  ffactionals  are  at  (0.25,  0.25)  and  (0.75,  0.75),  respectively.  The 
equilibrium  A-B  bond  distance  is  1.35  A,  and  the  angles  0  of  the  bonds  of  the  two  A-B  molecules  relative 
to  crystal  x-axis  are  ±29.1°. 

3.  DETAILS  OF  THE  CALCULATIONS 

3. 1  Cell-Linked  Lists.  The  two  types  of  simulations  that  we  report  here  require  calculating  the  energy 
and  energy  first  derivatives  of  Equation  (6).  In  principal,  Equation  (6)  requires  that  (N-l)  interactions 
must  be  calculated  for  each  atom  in  the  N-atom  system.  However,  for  a  specific  atom  i,  there  are  only 
a  small  number  of  neighbors  j  within  the  interaction  range  of  the  potential  and  therefore  are  the  only  ones 
out  of  the  N-atom  system  that  need  to  be  considered.  It  is  a  CPU-consuming  task  to  determine  by  brute 
force  which  of  the  (N-l)  atoms  are  within  the  interaction  range  of  atom  i.  This  problem  was  circumvented 
by  our  use  of  the  method  of  linked-lists,  described  in  detail  by  Allen  and  Tildesley  (1987).  This  method 
efficiently  sorts  the  atoms  into  indexed  bins  according  to  geometric  positioa  Each  bin  must  be  no  smaller 
than  the  cutoff  radius  of  the  interaction  potential.  The  scheme  ensures  that  for  a  specific  atom  in  a  bin, 
only  those  atoms  in  the  same  or  nearest-neighbor  bins  are  included  in  the  summation  of  Equation  (6).  The 
distance  between  atom  i  and  any  other  atom  not  in  the  same  or  nearest-neighbor  bins  (and  therefore, 
outside  of  the  interaction  range)  is  not  calculated,  thus  providing  considerable  CPU  savings.  For  the 
Monte  Carlo  simulations  (section  3.2),  each  bin  consists  of  nine  unit  cells  (three  in  each  of  the  x  and 
y  directions).  For  the  molecular  dynamics  simulations  that  are  used  in  calculating  the  equation  of  state 
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Table  1.  Parameters  and  Functional  Fonns  Used  for  the  Potential  Energy  Expression  in  Equation  (6) 
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molecular  bond  relative  to  the  crystal  x-axis. 


(section  3.3),  for  v/v0  greater  than  0.70,  each  bin  consists  of  four  unit  cells  (two  in  each  of  the  x  and 
y  directions).  For  v/v0  less  than  0.70,  each  bin  consists  of  six  unit  cells  (three  in  the  x  direction  and  two 
in  the  y  direction).  For  the  molecular  dynamics  computer  experiments  in  which  flyer  plate  impact  initiates 
the  shock  wave  (section  3.4),  each  bin  consists  of  four  unit  cells  (two  in  each  of  the  x  and  y  directions). 

3.2  NPT  Monte  Carlo  Simulations.  A  series  of  NPT  Monte  Carlo  simulations  were  performed  to 
determine  the  low-temperature  (20  K),  ambient  pressure  crystal  structure  and  sound  speed.  We  also 
performed  Monte  Carlo  calculations  for  pressures  up  to  0.55  eV/A2.  The  results  of  these  calculations  were 
used  as  starting  geometries  for  our  equation-of-state  calculations  detailed  as  follows  and  given  in  Table  2. 

The  sound  speed  of  a  crystal  is  proportional  to  the  slope  of  the  p-p  curve  (Fickett  and  Davis  1979; 
Fickett  1985);  to  obtain  the  sound  speed  of  our  low-temperature,  ambient  pressure  crystal,  we  calculated 
the  density  at  each  pressure  in  Table  2  and  fitted  the  pressure  from  0.015  to  0.0  eV/A2  to  a  cubic 
polynomial  in  density.  The  sound  speed  obtained  from  the  fit  is  1.2  km/s  for  this  crystal.  We  also  fitted 
<a>  and  <b>,  the  average  lattice  constants  of  the  unit  cell  in  the  x  and  y  directions,  respectively,  to 
quadratics  in  pressure,  to  provide  initial  lattice  parameters  for  the  trajectories  calculated  to  determine  the 
equation  of  state. 

The  simulation  box,  with  periodic  boundary  conditions  imposed  for  both  dimensions,  consisted  of  nine 
unit  cells  in  each  of  the  x  and  y  directions.  NPT  Monte  Carlo  simulations  were  symmetry-restricted:  In 
other  words,  at  each  attempted  step,  the  atomic  positions  of  all  atoms  in  the  system  were  determined  from 
the  positions  of  the  atoms  of  a  single  molecule,  which  we  will  denote  as  the  target  The  geometric 
parameters  sampled  in  these  simulations  were  the  lattice  parameters,  which  determine  the  volume  of  the 
unit  cell,  and  the  molecular  parameters  of  the  target  in  the  unit  cell.  The  position  of  the  COB  of  the  target 
molecule  was  fixed  at  the  lattice  fractional  (0.25,  0.25),  and  the  bond  length  and  orientation  angle  of  the 
A-B  molecule  relative  to  the  crystal  x-axis  was  allowed  to  be  varied  through  the  Monte  Carlo  sampling. 
Images  of  the  target  had  the  same  molecular  geometries  as  the  target,  but  were  translated  by  the  lattice 
spacings.  The  COB  of  the  second  molecule  in  the  unit  cell  as  the  target  was  fixed  at  the  lattice  fractional 
(0.75, 0.75),  and  the  molecule  was  assigned  the  same  bond  length  as  the  target,  but  had  the  negative  value 
of  the  orientation  angle  of  the  target  sampled  through  Monte  Carlo.  This  ensured  the  herringbone  lattice 
structure. 
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Table  2.  Lattice  Parameters,  Density,  and  Molecular  Geometry  vs.  Pressure 


Pressure 

<a> 

<b> 

<p> 

<rij> 

<eij> 

(eV/A2) 

(A) 

(A) 

(amu/A2) 

(A) 

O 

0.0 

4.34 

6.27 

4.4836 

1.349 

29.1 

0.00005 

4.33 

6.28 

4.4869 

1.349 

29.3 

0.0001 

4.33 

6.27 

4.4936 

1.349 

29.1 

0.0005 

4.32 

6.25 

4.5185 

1.349 

29.1 

0.001 

4.31 

6.22 

4.5505 

1.349 

29.0 

0.0025 

4.24 

6.19 

4.6476 

1.349 

29.8 

0.0125 

4.13 

5.94 

4.9735 

1.347 

29.0 

0.015 

4.11 

5.90 

5.0309 

1.346 

29.0 

0.05 

3.98 

5.67 

5.4054 

1.340 

29.3 

0.10 

3.85 

5.46 

5.8040 

1.332 

28.6 

0.15 

3.76 

5.33 

6.0878 

1.326 

29.1 

0.20 

3.71 

5.23 

6.2887 

1.321 

28.8 

0.25 

3.67 

5.14 

6.4687 

1.317 

28.6 

0.30 

3.63 

5.08 

6.6161 

1.314 

28.6 

0.35 

3.59 

5.02 

6.7703 

1.311 

28.5 

0.40 

3.56 

4.96 

6.9083 

1.310 

28.3 

0.45 

3.53 

4.91 

7.0398 

1.310 

28.1 

0.50 

3.50 

4.85 

7.1849 

1.312 

27.9 

0.55 

3.46 

4.77 

7.3939 

1.317 

27.5 

Several  Monte  Carlo  calculations  were  perforated  using  different  initial  configurations  to  see  if  the 
results  converged  to  the  same  value,  regardless  of  initial  state.  The  molecular  geometry  and  lattice 
spacings  for  the  crystal  were  initially  set  to  values  that  were  far  (up  to  ±33%)  from  equilibrium  values. 
For  example,  the  initial  bond  length  of  the  target  ranged  from  1.0  to  1.8  A,  and  the  lattice  spacings  were 
either  0.8  A  larger  or  smaller  than  equilibrium  values.  A  total  of  1,000  Markov  steps  were  attempted, 
during  which  time  the  system  relaxed  to  near  its  thermal  equilibrium  configuration.  At  this  point,  4,000 
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points  were  used  in  the  averaging  of  the  results.  Once  the  lattice  parameters  were  obtained,  we  determined 
the  equilibrium  orientation  of  the  molecules  in  the  crystal  corresponding  to  those  lattice  parameters  using 
the  Newton-Raphson  energy  minimization  method  (Press  et  al.  1986). 

3.3  Molecular  Dynamics  Simulations:  Equation-of-State  Calculations.  In  order  to  calculate  the 
Hugoniot  function  in  Equation  (5),  we  are  required  to  calculate  the  equation  of  state  of  the  system  behind 
the  shock  front  Additionally,  we  need  to  know  the  state  of  the  quiescent  system.  We  have  used 
molecular  dynamics  to  determine  these  states,  in  the  manner  outlined  by  Erpenbeck  (1992).  In  that  study, 
Erpenbeck  determined  the  equation  of  state  and  Hugoniot  curve  for  a  simple  diatomic  fluid  using 
thermodynamic  averages  calculated  from  molecular  dynamics  trajectory  ensembles.  The  initial  conditions 
for  each  trajectory  in  the  ensembles  were  selected  using  standard  Metropolis  Monte  Carlo  sampling.  We 
have  followed  the  procedure  outlined  in  Erpenbeck  (1992),  except  rather  than  use  ensembles  of  short-lived 
trajectories  for  our  averaging,  we  extracted  time-averaged  thermodynamic  properties  from  a  single 
long-lived  trajectory  for  a  given  set  of  initial  conditions.  Each  trajectory  was  integrated  until  the  averages 
of  the  thermodynamic  properties  converged.  Most  averages  converged  within  5  ps;  some  trajectories  were 
integrated  up  to  10  ps  before  convergence  was  met.  Thermodynamic  properties  were  calculated  for  v 
ranging  from  0.123  to  0.223  A2/amu.  We  also  differed  in  the  way  we  calculated  the  pressure  of  the 
system;  Erpenbeck  used  the  virial  function  to  calculate  his  pressures  (Erpenbeck  1992);  we  followed  the 
method  outlined  in  Tsai  (1979).  Periodic  boundary  conditions  were  imposed  in  both  x  and  y  directions. 
The  initial  conditions  of  the  crystal  for  a  simulation  corresponding  to  each  v  were  determined  from  the 
quadratic  fits  of  the  lattice  parameters  given  in  Table  2.  The  molecular  bond  length  and  orientation  of 
the  molecular  bond  relative  to  the  crystal  x-axis  were  set  to  1.35  A  and  29°,  respectively.  Kinetic  energy 
ranging  from  2,000  K  to  15,000  K  was  equipartitioned  among  the  atomic  momenta  components.  The 
equations  of  motion  of  the  system  were  then  integrated  for  approximately  6.5  ps.  The  final  conditions 
of  the  warmup  trajectory  were  then  used  as  the  initial  condition  for  the  trajectory  from  which  the 
thermodynamic  averages  would  be  extracted.  This  final  trajectory  was  integrated  until  the  thermodynamic 
averages  converged. 

Hamilton’s  equations  of  motion  for  this  system  were  integrated  using  an  Adams-Moulton  fourth-order 
predictor-corrector  (Miller  and  George  1972)  integrator  with  error  tolerance  set  to  10-5.  We  found  that 
the  results  using  this  tolerance  did  not  deviate  from  those  in  which  the  tolerance  was  set  much  smaller. 
Energy  conservation  was  monitored,  and  accuracy  to  0.0003  eV  was  obtained. 
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3.4  Molecular  Dynamics  Simulations:  Computer  Experiment  of  a  Detonatioa  We  utilized  a  method 
developed  by  Tsai  and  Trevino  (1981)  in  which  the  simulation  box  expands  into  the  undisturbed  region 
as  the  shock  wave  passes  throughout  the  crystal.  We  are  interested  only  in  the  region  immediately 
preceding  and  following  the  shock  front.  To  simulate  an  infinitely  large  crystal  several  microns  from  the 
shock  front  would  merely  increase  the  simulation  time,  without  increasing  our  knowledge  of  the 
phenomenon  of  interest,  namely  the  detonation  and  the  region  affected  by  it.  Therefore,  we  implemented 
the  following  scheme:  The  simulation  box  initially  consists  of  A-B  molecules  at  the  equilibrium  position. 
It  is  a  16  x  8  area  of  unit  cells,  with  periodic  boundary  conditions  imposed  in  the  y  direction  only. 
Figure  1  illustrates  the  initial  state  of  the  simulation  system.  For  discussion  purposes  throughout  this 
report,  we  will  denote  a  fragment  of  the  molecular  crystal  consisting  of  2  x  8  unit  cells  as  a  "slab."  To 
minimize  surface  effects  at  the  far  right  edge  of  the  cell  (which  is  furtherest  from  shock  impact  of  the 
plate),  we  held  a  slab  of  A-B  molecules  rigid  throughout  the  simulation,  with  the  molecules  fixed  in  their 
equilibrium  orientation.  All  other  A-B  molecules  and  flyer  plate  atoms  are  allowed  to  move  according 
to  the  equations  of  motion.  The  flyer  plate  is  a  slab  of  A-A  molecules  (located  at  the  far  left  of  Figure  1), 
chosen  because  of  their  stability  in  order  that  chemical  reactivity  would  not  contribute  to  the  mechanical 
energy  that  is  transferred  to  the  stationary  A-B  molecular  crystal  upon  impact.  Note  that  a  slab  in  this 
figure  is  highlighted  and  designated  as  "test  slab."  The  average  kinetic  energy  for  the  atoms  in  the  test 
slab  is  calculated  at  each  integration  step.  If  this  value  exceeds  15  K,  then  a  new  slab  of  A-B  molecules 
is  inserted  between  the  rigid  slab  and  the  far  right  edge  of  the  slab  of  atoms  that  are  allowed  to  move 
throughout  the  simulation.  The  molecules  in  the  new  slab  are  initially  at  their  equilibrium  position,  and 
kinetic  energy  corresponding  to  20  K  is  partitioned  between  the  x-  and  y-momentum  components  for  each 
atom.  When  a  new  slab  of  molecules  is  added,  the  test  slab  is  shifted  by  one  slab  length  to  the  right  In 
this  scheme,  the  length  of  the  undisturbed  crystal  is  constant  and  consists  of  seven  slabs  of  quiescent  A-B 
molecules  (not  including  the  rigid  slab).  We  found  that  the  energy  rapidly  equilibrated  in  this  region  and 
was  partitioned  equally  into  potential  and  kinetic  energy  (average  kinetic  energy  in  this  region  is  10  K). 
This  latter  observation  serves  as  ad  hoc  justification  for  the  treatment  of  the  undisturbed  region. 

Initial  conditions  were  selected  as  follows:  All  atoms  in  the  simulation  box  are  at  the  equilibrium 
position;  the  A-A  molecules  in  the  flyer  plate  have  the  same  lattice  parameters  and  orientational  angles 
as  the  A-B  molecules;  the  only  difference  in  crystal  structure  between  the  flyer  plate  and  quiescent  A-B 
crystal  is  the  molecular  size  of  the  molecules.  Each  atom  is  given  kinetic  energy  totaling  20  K,  partitioned 
between  the  x-  and  y-momentum  components.  A  short  warmup  trajectory  is  integrated  for  0.7  ps,  to  allow 
randomization  of  the  energy  in  the  crystal.  Because  the  flyer-plate  slab  of  A-A  molecules  is  not  in  the 
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equilibrium  position  (it  has  smaller  molecular  bond  distances  than  molecules  in  the  A-B  crystal),  some 
heat  will  be  released  into  Ihe  system  during  this  warmup  trajectory  as  the  A-A  atoms  relax  toward  the 
equilibrium  position  for  those  lattice  parameters.  We  found  through  monitoring  the  average  kinetic  energy 
of  the  A-A  molecules  in  the  flyer  plate  and  of  the  A-B  molecules  in  the  adjacent  slab  during  the  warmup 
trajectory  that  the  amount  of  heat  released  is  insignificant.  The  average  kinetic  energy  of  the  molecules 
in  both  slabs  fluctuated  about  10  K  throughout  the  warmup  trajectory.  After  the  warmup  trajectory,  the 
flyer  plate  atoms  are  given  impact  velocities  in  the  positive  x-direction.  As  the  equations  of  motion  of 
the  system  are  integrated,  the  flyer  plate  atoms  strike  the  quiescent  A-B  molecular  crystal.  We  found  that 
for  this  size  of  flyer  plate  (one  slab),  the  minimum  velocity  of  the  flyer  plate  to  initiate  an  unsupported 
detonation  is  4.7  km/s.  Anything  below  this  velocity  caused  a  few  reactions  at  the  initial  impact,  but 
apparently  were  not  enough  to  sustain  the  detonation. 

Energy  conservation  was  monitored,  and  accuracy  to  0.0003  eV  was  obtained  until  a  new  slab  of 
atoms  was  added.  At  this  point,  there  is  a  discontinuity  in  the  energy  of  the  system  (more  equations  of 
motion  are  being  integrated  due  to  the  addition  of  atoms),  but  energy  is  conserved  until  another  slab  is 
added  during  the  simulation. 

The  mass  density,  kinetic  temperature,  and  two-dimensional  pressures  in  the  steady  region  of  the 
reactive  flow  in  our  simulations  were  calculated  in  a  manner  outlined  in  Tsai  and  Trevino  (1981).  Our 
simulations  differ  from  the  piston-driven  shock  wave  results  calculated  by  Tsai  and  Trevino  (1981)  in  that 
our  simulations  result  in  unsupported  detonations.  An  unsupported  detonation  has  a  following  flow  behind 
the  steady  reaction  zone  that  changes  with  time,  and  the  conservation  equations  [Equations  (1M5)]  cannot 
be  applied  to  this  region  (Fickett  and  Davis  1979;  Fickett  1985).  Therefore,  we  had  to  make  some 
approximations  in  calculating  kinetic  temperature  and  two-dimensional  pressures  in  the  rarefaction  zone. 
To  calculate  the  kinetic  temperature  associated  with  thermal  motion  for  regions  throughout  the  crystal,  we 
had  to  remove  the  kinetic  energy  associated  with  mass  flow  velocity  for  the  region.  For  areas  within  the 
steady  reaction  zone,  the  local  mass  flow  can  be  determined  from  Equation  (1).  For  areas  in  the 
rarefaction  zone,  we  approximated  the  local  flow  velocity  to  be  the  center  of  mass  velocity  of  all  of  the 
particles  within  this  area. 

Equation  (2)  can  be  employed  to  calculate  the  pressure  through  the  shock  front  for  the  steady  portion 
of  the  shock  wave.  In  the  rarefaction  region.  Equation  (2)  no  longer  applies.  Instead,  the  instantaneous 
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stress  is  obtained  by  the  method  outlined  in  Tsai  (1979),  which  uses  the  forces  and  momentum  flux  that 
intercept  lines  moving  at  the  local  flow  velocity. 

4.  RESULTS 

4.1  Eauation-of-State  Calculations.  We  have  calculated,  using  molecular  dynamics,  the 
thermodynamic  properties  of  the  quiescent  A-B  molecular  crystal  and  the  sample  at  specific  volumes 
from  0.123  to  0.188  A2/amu  and  for  temperatures  ranging  from  approximately  2,000  K  up  to  10,000  K. 
The  results  are  shown  in  Table  3.  Figures  2  and  3  show  the  temperature  dependence  of  both  the  Hugoniot 
function  [Equation  (5)]  and  pressure  for  various  specific  volumes.  The  symbols  represent  the  results  of 
the  molecular  dynamics  averages  of  these  properties,  and  the  curves  in  Figure  3  are  fits  of  the  pressure 
to  quadratic  polynomials  in  temperature.  The  quadratic  fits  of  Equation  (5)  to  temperature  were  not 
illustrated  in  Figure  2,  for  clarity  of  the  figure.  The  Hugoniot  temperatures  (TH),  the  temperature  at  which 
Equation  (5)  is  zero  for  each  specific  volume,  were  extrapolated  from  the  quadratic  fits  of  Equation  (5) 
to  temperature,  and  are  listed  in  Table  4.  The  corresponding  Hugoniot  pressure,  PH>  was  calculated  at  the 
TH  for  each  specific  volume  using  the  quadratic  fit  of  pressure  to  temperature,  and  are  also  listed  in 
Table  4.  These  values  can  be  used  in  Equation  (3)  to  determine  detonation  velocities,  D,  as  functions  of 
specific  volume.  The  detonation  velocities  corresponding  to  each  TH  and  PH  are  given  in  Table  4  and 
shown  as  a  function  of  specific  volume  in  Figure  4.  The  Chapman-Jouguet  velocity  is  the  minimum 
detonation  velocity  corresponding  to  the  set  of  Hugoniot  pressures  and  specific  volumes.  The  curve  shown 
in  Figure  4  is  a  fit  of  the  detonation  velocities  to  a  quadratic  polynomial  in  specific  volume.  The  position 
of  the  minimum  of  D,  calculated  using  the  quadratic  polynomial  fit  shown  in  Figure  4,  corresponds  to  a 
specific  volume  of  0.141  A2/amu  (or  density  of  7.09  amu/A2).  The  detonation  velocity  at  this  density  is 
0.717  AA.u.  (7.0  km/s). 

Figure  5  shows  the  Hugoniot  curve  and  the  Rayleigh  line  that  uses  the  Chapman-Jouguet  detonation 
velocity  determined  from  the  polynomial  fit  shown  in  Figure  4.  Indeed  the  curves  intersect  at  the  C-J 
point  The  C-J  pressure  is  0.86  eV/A2. 

4.2  Molecular  Dynamics  Simulation  of  Detonation.  Molecular  dynamics  simulations  of  a  flyer  plate 
impacting  a  quiescent  crystal  were  performed  with  six  impact  velocities  ranging  from  4.6  km/s  to  12  km/s. 
The  simulations  will  be  denoted  hereafter  by  the  impact  velocity  of  the  flyer  plate.  The  position  of  the 
shock  wave  as  a  function  of  time  for  each  of  these  six  simulations  are  illustrated  in  Figure  6.  The 
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Table  3.  Thermodynamic  Properties  vs.  Specific  Volume 


y 

T 

e 

P 

(A2/amu) 

(K) 

(eV/amu) 

(eV/A2) 

0.223048 

10 

-0.016945 

0.000790 

0.187541 

7,702 

-0.005269 

0.468962 

8,541 

0.000381 

0.483094 

9,116 

0.006033 

0.478626 

10,076 

0.011683 

0.484239 

0.167262 

6,671 

-0.003577 

0.623313 

7,302 

0.002074 

0.625804 

8,303 

0.007724 

0.618631 

8,845 

0.013375 

0.617379 

0.155948 

4,247 

-0.007119 

0.648500 

5,763 

-0.001468 

0.707511 

5,836 

-0.001468 

0.721500 

6,491 

0.004183 

0.705964 

6,607 

0.004183 

0.725363 

7,197 

0.009833 

0.703152 

7,444 

0.009832 

0.697245 

1 

0.015484 

0.015484 

0.681139 

0.716467 

0.147425 

4,215 

-0.004544 

0.756504 

5,026 

0.001107 

0.778300 

5,791 

0.006758 

0.753580 

6,923 

0.012409 

0.797745 

7,689 

0.018059 

0.786444 

0.138170 

3,699 

-0.001146 

0.871516 

4,645 

0.005070 

0.862671 

5,311 

0.010720 

0.846950 

6,531 

0.016370 

0.889304 

7,386 

0.022021 

0.879788 

0.127342 

2,608 

0.002651 

0.919353 

3,108 

0.005477 

0.922500 

3,962 

0.011127 

0.900155 

4,806 

0.016778 

0.906584 

5,816 

0.022429 

0.931453 

6,999 

0.028079 

0.983333 

0.123003 

2,876 

0.008239 

0.918458 

3,704 

0.013890 

0.941040 

4,661 

0.019540 

0.942793 

5,672 

0.025191 

0.971843 

6,927 

0.030842 

1.021582 
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Figure  3.  Pressure  vs.  temperature  for  various  specific  volumes,  values  of  which  are  denoted  in  the 
legend.  These  are  the  results  of  molecular  dynamics  equation-of-state  calculations.  The  solid 
lines  denote  Quadratic  fits  in  temperature. 
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Table  4.  Hugoniot  Temperatures,  Pressures,  and  Detonation  Velocities  vs.  Specific  Volume 


V 

th 

PH 

D 

(A2/amu) 

(K) 

(eV/A2) 

(A/t.u.)a 

0.187541 

7,335 

0.46463 

0.806168 

0.167262 

7,235 

0.62413 

0.745588 

0.155948 

7,063 

0.71394 

0.727154 

0.147425 

6,834 

0.78271 

0.717220 

0.138170 

7,125 

0.88206 

0.718711 

0.127342 

7,500 

1.01308 

0.725406 

0.123003 

8,051 

1.07319 

0.730263 

1  tu.  =  1.018066  x  1(T14  s. 


Figure  4.  Detonation  velocity  D  as  a  function  of  specific  volume  obtained  from  the  Hugoniot 
temperatures  and  pressures  at  each  specific  volume  (see  Equation  2).  The  solid  curve  denotes 
a  quadratic  fit  in  specific  volume.  One  t.u.  equals  1.018066  x  10  14  s.  ~~~ 
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eauation-of-state  calculations.  The  straight  line  tangent  to  the  Hugoniot  curve  is  the  Rayleigh 
line  corresponding  to  the  Chapman-Jouguet  condition.  The  slope  of  this  line  is  consistent  with 
the  Chapman-Jouguet  velocity  determined  from  the  data  shown  in  Figure  4. 


Figure  6.  Position  of  the  shock  front  as  a  function  of  time  for  the  six  molecular  dynamics  simulations. 
The  numbers  denote  the  velocities  of  the  flyer  plate  in  kilometers  per  seconds. 
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threshold  for  initiation  of  detonation  for  this  plate  thickness  is  found  to  be  4.7  km/s;  detonation  was  not 
sustained  for  plate  impacts  smaller  than  this.  With  the  exception  of  the  4.6  km/s  simulation,  the  slopes 
of  the  shock  fronts  are  the  same  by  3.5  ps  into  each  simulation.  Before  this  time,  the  slopes  differ  as  the 
detonations  reach  steady  state.  The  steady-state  detonation  velocity  is  6.6  km/s.  The  faster  the  impact 
velocity,  the  sooner  the  detonation  reaches  steady  state.  The  shock  wave  initiated  with  flyer  plate  impact 
of  12  km/s  reaches  the  steady-state  detonation  velocity  within  0.1  ps.  The  detonation  velocity  determined 
through  these  computer  experiments  is  6.1%  smaller  than  the  C-J  detonation  velocity  predicted  from 
hydrodynamic  theory  (see  section  4.1). 

Thermodynamic  property  profiles  for  the  computer  experiments  have  similar  features.  A  typical 
snapshot  of  the  system  for  the  12  km/s  simulation  at  7.8  ps  is  shown  in  Figure  7;  thermodynamic  property 
and  species  profiles  corresponding  to  this  time  are  shown  in  Figure  8.  Figure  7  shows  three  distinct 
regions:  Undisturbed  crystal,  the  reaction  zone,  in  which  the  molecules  are  compressed  and  are 
undergoing  reaction,  and  the  rarefaction  region,  which  consists  of  vibrationally-excited  homonuclear 
products.  We  have  defined  the  reaction  zone  as  the  area  between  the  position  of  the  shock  front  and  the 
point  (along  the  x-axis)  at  which  the  number  of  reacted  molecules  (dissociated  from  original  molecular 
partner)  exceeds  the  number  of  unreacted  A-B  molecules.  The  point  at  which  the  number  of  reactions 
exceeds  the  number  of  unreacted  atoms  is  illustrated  in  Figure  8(d);  it  is  approximately  14  A  behind  the 
shock  front  in  this  snapshot.  We  have  found  that  the  size  of  the  reaction  zone  is  steady  in  time,  and  is 
on  average  14  A  in  width.  This  is  illustrated  in  Figure  9,  which  shows  the  width  of  the  reaction  zone 
throughout  the  12  km/s  simulation.  Additionally,  the  composition  and  properties  of  this  zone  are  steady 
in  time.  The  thermodynamic  properties  in  this  region  are  those  we  wish  to  compare  with  the 
hydrodynamic  predictions.  We  have  calculated  the  thermodynamic  properties  of  a  thin  area  of  the  sample 
directly  behind  the  shock  front  over  the  life  of  the  12  km/s  trajectory.  This  area  has  the  same  dimensions 
of  a  "slab"  as  defined  in  section  3.4.  Figures  10(a-c)  show  the  density,  pressure,  and  kinetic  temperature 
of  this  region  over  time.  The  properties  are  well  behaved  and  steady.  The  average  temperature,  pressure 
and  density  of  this  region  are:  7,435  K,  0.88  eV/A2,  and  7.13  amu/A2.  The  average  pressure  and  specific 
volume  of  these  regions  differ  by  2.6%  and  0.6%,  respectively,  from  the  hydrodynamic  predictions  of  the 
C-J  values.  The  agreement  between  the  two  theories  is  good. 

The  differences  in  the  results  of  the  two  calculations,  although  small,  can  be  attributed  to  differences 
in  state  composition.  In  the  equation-of-state  calculations,  an  equimolar  mixture  of  A  and  B  atoms  was 
used.  In  the  detonation  simulations,  however,  the  number  density  of  B  atoms  in  the  slab  behind  the  shock 
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snapshot  are  300  x  50  A. 


Figure  8.  Thermodynamic  property  profiles  of  the  system  corresponding  to  the  conditions  described  in 
Figure  7.  Solid  and  dashed  curves  in  (d)  denote  unreacted  and  reacted  A-B  molecules. 
respectively.  The  smooth  curves  in  (aHc)  are  guides  to  the  eve. 
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Figure  9.  Reaction  zone  length  as  a  function  of  time  through  the  12  km/s  simulation. 

front  was  consistently  higher  than  the  number  density  of  A  atoms.  The  mixture  behind  the  shock  front 
consists,  on  average,  of  57%  B  atoms,  43%  A  atoms.  Therefore,  although  the  C-J  density  agrees  well  with 
the  density  behind  the  shock  front  in  the  computer  experiment,  the  chemical  composition  between  the  two 
situations  differs.  The  larger  percentage  of  the  heavier  B  atoms  in  the  detonation  simulation  is  consistent 
with  a  detonation  velocity  lower  than  the  C-J  detonation  velocity  predicted  from  the  results  of  the 
equation-of-state  calculations. 

The  zone  behind  the  reaction  zone  (the  rarefaction  zone)  is  unsteady;  its  properties  change  with  time. 
The  configuration  of  the  system  produced  in  our  computer  experiments  appears  to  be  representative  of  an 
unsupported  detonation  as  described  in  Fickett  and  Davis  (1979)  and  Fickett  (1985),  and  is  almost 
described  by  the  simplest  theory  (Fickett  and  Davis  1979;  Fickett  1985).  The  flow  appears  to  be 
one-dimensional,  the  detonation  front  is  almost  a  jump  discontinuity,  and  the  material  emerging  from  the 
front  (the  reaction  zone)  is  independent  of  time.  Additionally,  the  model  used  in  these  computer 
experiments  provides  almost  instantaneous  reaction,  with  only  slightly  more  than  14  A  between  the  shock 
front  and  complete  reaction.  These  properties,  as  well  as  the  good  agreement  with  the  hydrodynamic 
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Figure  10.  Densities,  pressures,  and  temperatures  as  functions  of  time  for  simulations  with  flyer-plate  velocities  of  12  km/s  (aMc).  4.7  km/s 
(d)-(e).  and  4.6  km/s  (g)-ffl.  Note  the  significant  differences  in  the  temperature  scales. 


predictions,  indicate  that  both  the  method  of  molecular  dynamics  and  our  model  of  an  energetic  crystal 
can  reasonably  describe  the  phenomenon  of  detonation. 

We  have  also  calculated  the  time  dependence  of  the  thermodynamic  properties  of  a  slab  of  material 
(as  defined  in  section  3.4)  immediately  behind  the  shock  fronts  in  the  computer  experiments  for  plate 
impacts  of  4.7  and  4.6  km/s.  The  thermodynamic  properties  behind  the  shock  fronts  for  these  two 
simulations  were  evaluated  as  for  the  12  km/s  simulation  [see  Figures  10(d-i)].  The  4.6  km/s  system 
reaches  a  maximum  kinetic  temperature,  pressure  and  density  of  3,800  K,  0.61  eV/A2  and  7.6  amu/A2, 
respectively,  within  1  ps  into  the  simulation,  but  these  values  are  not  maintained  and  rapidly  drop  to 
significantly  lower  values  as  the  simulation  progresses.  The  maximum  pressure  reached  in  this  experiment 
is  well  below  the  Chapman- Jouguet  pressure,  although  the  C-J  density  is  reached  (but  not  maintained) 
during  this  trajectory.  The  profiles  for  the  4.7  km/s  simulation  are  very  similar  to  those  of  the  4.6  km/s 
simulation  for  the  first  two  picoseconds;  however,  the  theimodynamic  properties  approach  the  steady-state 
detonation  averages  after  this  time. 

It  is  worthwhile  to  examine  more  closely  the  differences  in  the  profiles  between  the  simulations  of 
the  4.6  km/s  and  4.7  km/s  flyer-plate  simulations  since  the  former  does  not  lead  to  sustained  detonation 
and  the  latter  does.  Figure  11  shows  the  densities,  pressures,  and  temperatures  of  the  reaction  zone  as 
functions  of  time  during  the  first  two  picoseconds  of  both  the  4.6  km/s  and  4.7  km/s  simulations.  The 
main  differences  in  these  properties  for  the  two  simulations  occur  between  0.6-1 .0  ps,  and  then  after 
1.5  ps.  At  0.7  ps,  there  is  a  sharp  increase  in  the  kinetic  temperature  for  the  4.7  km/s  simulation  that  does 
not  occur  for  the  4.6  km/s  simulation.  The  temperature  then  fluctuates  near  3,000  K  for  0.3  ps  in  the 
4.7  km/s  simulation,  whereas  the  temperature  for  the  4.6  km/s  simulation  falls  below  2,000  K  during  this 
same  time.  The  pressure  in  the  4.7  km/s  simulation  averages  approximately  0.55  eV/A2  during  this  time, 
whereas  the  pressure  of  the  4.6  km/s  simulation  falls  to  half  that  value.  Finally,  the  density  of  the 
4.7  km/s  simulation  from  0.6  ps  to  1.0  ps  fluctuates  near  the  C-J  value,  but  the  density  for  the  4.6  km/s 
simulation  during  this  time  period  is  well  below  the  C-J  value. 

The  two  simulations  predict  similar  behavior  in  properties  from  1.0  ps  to  1.5  ps,  but  then  the  curves 
diverge.  All  curves  corresponding  to  the  4.6  km/s  simulation  decrease  monotonically.  For  the  4.7  km/s 
simulation,  the  C-J  density  is  reached  and  subsequently  maintained  at  1.6  ps,  well  before  the  C-J  pressure 
is  reached  0.6  ps  later.  The  temperature  for  the  4.7  km/s  simulation  does  not  reach  a  steady  value  until 
6  ps  into  the  trajectory;  it  monotonically  increases  to  5,500  K  at  approximately  2.2  ps,  and  fluctuates  about 
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Figure  11.  A  comparison  of  (a)  density,  (b)  pressure,  and  (c)  temperature  as  functions  of  time  for  the  two 
conditions  of  flyer-plate  velocity  which  bracket  the  threshold  for  sustained  unsupported 


detonation.  The  straight  lines  in  the  density  and  pressure  figures  denote  the  C-J  values  for 


these  quantities. 
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this  value  for  approximately  2  ps.  It  then  increases  to  6,500  K  for  another  2  ps,  at  which  point  it  appears 
to  increase  to  7,200  K,  about  which  it  fluctuates  for  the  remainder  of  the  trajectoiy.  These  fluctuations 
throughout  the  simulations,  particularly  at  the  beginning,  are  suggestive  of  a  system  sampling  a  critical 
region  of  phase  space  separating  two  distinctly  different  results:  namely,  sustained,  unsupported 
detonation  vs.  nonreactive  shock.  From  these  results,  it  would  appear  that  attainment  of  the  C-J  density 
from  shock  impact  is  the  determining  factor  for  the  sustenance  of  a  detonation. 

5.  CONCLUSIONS 

We  have  presented  a  comparative  study  of  molecular  dynamics  computer  experiments  of  an 
unsupported  detonation  to  hydrodynamic  predictions  based  on  the  classical  conservation  equations  that 
relate  mass,  momentum,  and  energy  of  the  quiescent  crystal  with  the  state  behind  the  detonation  wave 
(Fickett  and  Davis  1979;  Fickett  1985).  We  calculated  the  equation  of  state  of  the  system,  through 
molecular  dynamics,  in  order  to  evaluate  the  classical  conservation  equations  and  generate  the  Hugoniot 
curve  for  this  system.  The  model  used  in  both  the  equation-of-state  calculations  and  in  the  computer 
experiments  describes  a  reactive  crystal  consisting  of  heteronuclear  diatomic  molecules  that  release  heat 
upon  formation  of  the  homonuclear  diatomic  products.  All  calculations  presented  herein  are 
two-dimensional. 

The  equation  of  state  of  an  equimolar  mixture  of  A  and  B  atoms  is  determined  from  thermodynamic 
averages  obtained  through  molecular  dynamics  simulations  for  specific  volumes  ranging  from  v0,  the 
low-pressure  reduced  volume,  to  0.55  v0.  Equilibrium  temperatures  and  pressures  were  determined  for 
each  reduced  volume,  and  the  Hugoniot  curve  was  produced.  The  Chapman-Jouguet  state  was  then 
determined;  the  C-J  detonation  velocity,  density,  and  pressure  are  predicted  to  be  7.0  km/s,  7.09  amu/A2, 
and  0.86  eV/A2,  respectively. 

The  computer  experiments  simulate  shock-initiated  (through  flyer-plate  impact)  reactions  in  a  model 
energetic  crystal.  For  plate  impacts  with  velocities  no  less  than  4.7  km/s,  the  shock  front  and  reaction 
zone  propagate  through  the  crystal  at  a  steady  rate  of  6.6  km/s.  The  thermodynamic  properties  and  width 
of  the  reaction  zone  are  time-independent  The  average  width  of  the  reaction  zone  is  14  A;  this  width  is 
the  same  for  all  plate  impact  velocities  no  less  than  4.7  km/s.  Time-averaged  density,  pressure,  and 
temperature  immediately  behind  the  shock  front  in  the  reaction  zone  are  7.13  amu/A2,  0.88  eV/A2,  and 
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7,435  K,  respectively.  The  following  flow  is  time-dependent;  the  properties  and  width  of  this  region 
change  as  the  simulations  progress  in  time. 

Agreement  between  the  computer  experiment  and  the  hydrodynamic  predictions  is  good.  The  latgest 
discrepancy  is  a  6%  difference  in  the  detonation  velocities,  which  we  attribute  to  differences  in  the 
chemical  composition  of  the  system  used  in  the  equation-of-state  calculation  and  that  of  the  reaction  zone 
behind  the  shock  front  in  the  computer  experiments.  The  equation-of-state  calculations  had  an  equimolar 
distribution  of  A  and  B  atoms  in  the  system,  but  the  state  behind  the  shock  front  in  the  computer 
experiment  had,  on  average,  a  larger  number  density  of  B  atoms  (57%  of  the  total  number)  than  A  atoms 
(43%  of  the  total  number)  over  the  lifetime  of  the  simulation.  The  larger  concentration  of  the  heavier  B 
atoms  in  the  computer  experiment  could  explain  the  smaller  detonation  velocity  observed  in  the  computer 
experiments.  The  C-J  pressure  and  density  are  in  remarkably  good  agreement  with  the  computer 
simulations  (within  2.5%)  even  though  the  chemical  composition  behind  the  shock  front  in  the  computer 
experiment  differs  from  that  used  in  the  equation-of-state  calculations. 

Thermodynamic  properties  of  thin  regions  immediately  behind  the  shock  fronts  were  monitored  in  time 
for  computer  experiments  in  which  flyer  plates  strike  the  quiescent  molecular  crystal  with  velocities  of 
4.6  km/s  and  4.7  km/s,  respectively.  Detonation  is  not  sustained  for  the  4.6  km/s  impact;  an  unsupported 
detonation  results  from  the  4.7  km/s  impact.  For  the  first  1.5  ps  of  both  simulations,  the  thermodynamic 
properties  have  similar  values.  After  this  time,  the  behavior  of  the  properties  of  the  two  simulations 
diverge.  All  properties  corresponding  to  the  simulation  with  4.6  km/s  flyer-plate  impact  decrease 
monotonically  after  1.5  ps.  At  1.6  ps,  the  density  of  the  system  corresponding  to  the  4.7  km/s  flyer-plate 
impact  reaches  and  maintains  the  C-J  value.  The  pressure  for  the  4.7  km/s  simulation  subsequently 
increases  monotonically  to  the  C-J  value  only  after  the  C-J  density  is  attained,  suggesting  that  if  the  C-J 
density  is  reached,  then  a  detonation  will  be  sustained. 

These  results  show  that  the  method  of  molecular  dynamics  and  our  model  of  a  reactive  energetic 
molecular  crystal  can  be  used  to  reasonably  simulate  the  phenomenon  of  detonation. 
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ATTN  SARD  TT  MR  J  APPEL 
PENTAGON 

WASHINGTON  DC  20310-0103 

1  HQDA  OASA  RDA 

ATTN  DR  C  H  CHURCH 
PENTAGON  ROOM  3E486 
WASHINGTON  DC  20310-0103 

4  COMMANDER 

US  ARMY  RESEARCH  OFFICE 
ATTN  R  GHIRARDELLI 
D  MANN 
R  SINGLETON 
R  SHAW 
P  O  BOX  12211 

RSCH  TRNGLE  PK  NC  27709-2211 

1  DIRECTOR 

ARMY  RESEARCH  OFFICE 
ATTN  AMXRO  RT  IP  LIB  SERVICES 
P  O  BOX  12211 

RSCH  TRNGLE  PK  NC  27709-2211 

2  COMMANDER 
US  ARMY  ARDEC 

ATTN  SMCAR  AEE  B  D  S  DOWNS 
PCTNY  ARSNL  NJ  07806-5000 

2  COMMANDER 

US  ARMY  ARDEC 

ATTN  SMCAR  AEE  J  A  LANNON 

PCTNY  ARSNL  NJ  07806-5000 

1  COMMANDER 
US  ARMY  ARDEC 

ATTN  SMCAR  AEE  BR  L  HARRIS 
PCTNY  ARSNL  NJ  07806-5000 

2  COMMANDER 

US  ARMY  MISSILE  COMMAND 
ATTN  AMSMI  RD  PR  E  A  R  MAYKUT 
AMSMI  RD  PR  P  R  BETTS 
REDSTONE  ARSENAL  AL  35809 


1  OFFICE  OF  NAVAL  RESEARCH 
DEPARTMENT  OF  THE  NAVY 
ATTN  R  S  MILLER  CODE  432 
800  N  QUINCY  STREET 
ARLINGTON  VA  22217 

1  COMMANDER 

NAVAL  AIR  SYSTEMS  COMMAND 
ATTN  J  RAMNARACE  AIR  5411 1C 
WASHINGTON  DC  20360 

2  COMMANDER 

NAVAL  SURFACE  WARFARE  CENTER 
ATTN  R  BERNECKER  R  13 
G  B  WILMOT  R  16 
SILVER  SPRING  MD  20903-5000 

5  COMMANDER 

NAVAL  RESEARCH  LABORATORY 
ATIN  M  C  LIN 
J  MCDONALD 
E  ORAN 
JSHNUR 

R  J  DOYLE  CODE  6110 
WASHINGTON  DC  20375 

2  COMMANDER 

NAVAL  WEAPONS  CENTER 
ATTN  T  BOGGS  CODE  388 
T  PARR  CODE  3895 
CHINA  LAKE  CA  93555-6001 

1  SUPERINTENDENT 

NAVAL  POSTGRADUATE  SCHOOL 
DEPT  OF  AERONAUTICS 
ATTN  D  W  NETZER 
MONTEREY  CA  93940 

3  ALLSCF 
ATTN  R  CORLEY 
R  GEISLER 

J  LEVINE 

EDWARDS  AFB  CA  93523-5000 

1  AFOSR 

ATTN  J  M  TISHKOFF 
BOLLING  AIR  FORCE  BASE 
WASHINGTON  DC  20332 
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1  OSD  SDIO  1ST 

ATTN  L  CAVENY 
PENTAGON 

WASHINGTON  DC  20301-7100 

1  COMMANDANT 

USAFAS 

ATTN  ATSF  TSM  CN 
FORT  SILL  OK  73503-5600 

1  UNTV  OF  DAYTON  RSCH  INSTITUTE 

ATTN  D  CAMPBELL 
ALPAP 

EDWARDS  AFB  CA  93523 

1  NASA 

LANGLEY  RESEARCH  CENTER 
ATTN  G  B  NORTH  AM  MS  168 
LANGLEY  STATION 
HAMPTON  VA  23365 

4  NATIONAL  BUREAU  OF  STANDARDS 
US  DEPARTMENT  OF  COMMERCE 
ATTN  J  HASTE 
M JACOX 
T  KASHIWAGI 
..  H  SEMERJIAN 

WASHINGTON  DC  20234 

2  DIRECTOR 

LAWRENCE  LIVERMORE  NATIONAL  LAB 
ATTN  C  WESTBROOK 
W  TAO  MS  L  282 
P  O  BOX  808 
LIVERMORE  CA  94550 

1  DIRECTOR 

LOS  ALAMOS  NATIONAL  LAB 
ATTN  B  NICHOLS  T7  MS  B284 
P  O  BOX  1663 
LOS  ALAMOS  NM  87545 

2  PRINCETON  COMBUSTION 
RESEARCH  LABORATORES  INC 
ATTN  N  A  MESSINA 

M  SUMMERFELD 
PRINCETON  CORPORATE  PLAZA 
BLDG  IV  SUITE  119 
11  DEERPARK  DRIVE 
MONMOUTH  JUNCTION  NJ  08852 


3  DIRECTOR 

SANDIA  NATIONAL  LABORATORES 
DIVISION  8354 
ATTN  S  JOHNSTON 
P  MATTERN 
D  STEPHENSON 
LIVERMORE  CA  94550 

1  BRIGHAM  YOUNG  UNIVERSITY 

DEPT  OF  CHEMICAL  ENGINEERING 
ATTN  M  W  BECKSTEAD 
PROVO  UT  84058 

1  CALEORNIA  INSTITUTE  OF  TECH 

JET  PROPULSION  LABORATORY 
ATTN  L  STRAND  MS  125  224 
4800  OAK  GROVE  DRIVE 
PASADENA  CA  91109 

1  CALEORNIA  INSTITUTE  OF  TECHNOLOGY 

ATTN  F  E  C  CULICK  MC  301  46 
204  KARMAN  LAB 
PASADENA  CA  91125 

1  UNIVERSITY  OF  CALEORNIA 

LOS  ALAMOS  SCIENTIFIC  LAB 
P  O  BOX  1663  MAIL  STOP  B216 
LOS  ALAMOS  NM  87545 

1  UNIVERSITY  OF  CALEORNIA  BERKELEY 

CHEMISTRY  DEPARMENT 
ATTN  C  BRADLEY  MOORE 
211  LEWIS  HALL 
BERKELEY  CA  94720 

1  UNIVERSITY  OF  CALEORNIA  SAN  DIEGO 
ATTN  F  A  WILLIAMS 

AMES  B010 
LA  JOLLA  CA  92093 

2  UNIV  OF  CALEORNIA  SANTA  BARBARA 
QUANTUM  INSTITUTE 

ATTN  K  SCHOFIELD 

M  STEINBERG 

SANTA  BARBARA  CA  93106 

1  UNIV  OF  COLORADO  AT  BOULDER 

ENGINEERING  CENTER 
ATTN  J  DAILY 
CAMPUS  BOX  427 
BOULDER  CO  80309-0427 
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3  UNIV  OF  SOUTHERN  CALIFORNIA 

DEPT  OF  CHEMISTRY 
ATTN  R  BEAUDET 
S  BENSON 
C  WTTTIG 

LOS  ANGELES  CA  90007 

1  CORNELL  UNIVERSITY 

DEPARTMENT  OF  CHEMISTRY 
ATTN  T  A  COOL 
BAKER  LABORATORY 
ITHACA  NY  14853 

1  UNIVERSITY  OF  DELAWARE 

CHEMISTRY  DEPARTMENT 
ATTN  T  BRILL 
NEWARK  DE  19711 

1  UNIVERSITY  OF  FLORIDA 

DEPT  OF  CHEMISTRY 
ATTN  J  WINEFORDNER 
GAINESVILLE  FL  32611 

3  GEORGIA  INSTITUTE  OF  TECHNOLOGY 

SCHOOL  OF  AEROSPACE  ENGINEERING 
ATTN  E  PRICE 
W  C  STRAHLE 
BTZEMN 

ATLANTA  GA  30332 

1  UNIVERSITY  OF  ILLINOIS 

DEPT  OF  MECH  ENG 
ATTN  H  KRIER 
144MEB  1206  W  GREEN  ST 
URBANA  IL  61801 

1  THE  JOHNS  HOPKINS  UNIV  CPIA 

ATTN  T  W  CHRISTIAN 
10630  LITTLE  PATUXENT  PKWY 
SUITE  202 

COLUMBIA  MD  21044-3200 

1  UNIVERSITY  OF  MICHIGAN 

GAS  DYNAMICS  LAB 
ATTN  G  M  FAETH 
AEROSPACE  ENGINEERING  BLDG 
ANN  ARBOR  MI  48109-2140 

1  UNIVERSITY  OF  MINNESOTA 

DEPT  OF  MECHANICAL  ENGINEERING 
ATTN  E  FLETCHER 
MINNEAPOLIS  MN  55455 


4  PENNSYLVANIA  STATE  UNIVERSITY 
DEPT  OF  MECHANICAL  ENGINEERING 
ATTN  K  KUO 
MMICCI 
S  THYNELL 
V  YANG 

UNIVERSITY  PARK  PA  16802 

2  PRINCETON  UNIVERSITY 

FORRESTAL  CAMPUS  LIBRARY 
ATTN  K  BREZINSKY 
I GLASSMAN 
P  O  BOX  710 
PRINCETON  NJ  08540 

1  PURDUE  UNIVERSITY 

SCHL  OF  AERONAUTICS  &  ASTRONAUTICS 

ATTN  J  R  OSBORN 

GRISSOM  HALL 

WEST  LAFAYETTE  IN  47906 

1  PURDUE  UNIVERSITY 
DEPARTMENT  OF  CHEMISTRY 
ATTN  E  GRANT 

WEST  LAFAYETTE  IN  47906 

2  PURDUE  UNIVERSITY 

SCHL  OF  MECHANICAL  ENGNRNG 
ATTN  N  M  LAURENDEAU 
SNBMURIHY 
TSPC  CHAFFEE  HALL 
WEST  LAFAYETTE  IN  47906 

1  RENSSELAER  POLYTECHNIC  INST 

DEPT  OF  CHEMICAL  ENGINEERING 
ATTN  A  FONTUN 
TROY  NY  12181 

1  STANFORD  UNIVERSITY 

DEPT  OF  MECHANICAL  ENGINEERING 
ATTN  R  HANSON 
STANFORD  CA  94305 

1  UNIVERSITY  OF  TEXAS 

DEPT  OF  CHEMISTRY 
ATTN  W  GARDINER 
AUSTIN  TX  78712 

1  VA  POLYTECH  INST  AND  STATE  UNIV 

ATTN  J  A  SCHETZ 
BLACKSBURG  VA  24061 
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APPLIED  COMBUSTION  TECHNOLOGY  INC 
ATTN  A  M  VARNEY 
P  O  BOX  607885 
ORLANDO  FL  32860 

APPLIED  MECHANICS  REVIEWS 
ASME 

ATTN  R  E  WHITE  &  A  B  WENZEL 
345  E  47TH  STREET 
NEW  YORK  NY  10017 

TEXTRON  DEFENSE  SYSTEMS 
ATTN  A  PATRICK 
2385  REVERE  BEACH  PARKWAY 
EVERETT  MA  02149-5900 

BATTELLE 

TWSTIAC 

505  KING  AVENUE 

COLUMBUS  OH  43201-2693 

COHEN  PROFESSIONAL  SERVICES 
ATTN  N  S  COHEN 
141  CHANNING  STREET 
REDLANDS  CA  92373 

EXXON  RESEARCH  &  ENG  CO 

ATTN  A  DEAN 

ROUTE  22E 

ANN  AND  ALE  NJ  08801 

GENERAL  APPLIED  SCIENCE  LABS  INC 
77  RAYNOR  AVENUE 
RONKONKAMA  NY  11779-6649 

GENERAL  ELECTRIC  ORDNANCE  SYSTEMS 
ATTN  JMANDZY 
100  PLASTICS  AVENUE 
PITTSFIELD  MA  01203 

GENERAL  MOTORS  RSCH  LABS 
PHYSICAL  CHEMISTRY  DEPARTMENT 
ATTN  T  SLOANE 
WARREN  MI  48090-9055 

HERCULES  INC 
ATTN  W  B  WALKUP 
E A  YOUNT 
P  O  BOX  210 

ROCKET  CENTER  WV  26726 


1  HERCULES  INC 

ATTN  R  V  CARTWRIGHT 
100  HOWARD  BLVD 
KENVIL  NJ  07847 

1  ALLIANT  TECHS YSTEMS  INC 

MARINE  SYSTEMS  GROUP 
ATTN  D  E  BRODEN  MS  MN50  2000 
600  2ND  STREET  NE 
HOPKINS  MN  55343 

1  ALLIANT  TECHS  YSTEMS  INC 

ATTN  R  E  TOMPKINS 
MN  11  2720 

600  SECOND  ST  NORTH 
HOPKINS  MN  55343 

1  IBM  CORPORATION 

RESEARCH  DIVISION 
ATTN  A  C  TAM 
5600  COTTLE  ROAD 
SAN  JOSE  CA  95193 

1  HT  RESEARCH  INSTITUTE 

ATTN  R  F  REMALY 
10  WEST  35TH  STREET 
CHICAGO  IL  60616 

1  LOCKHEED  MISSILES  &  SPACE  CO 

ATTN  GEORGE  LO 
3251  HANOVER  STREET 
DEPT  52  35  B204  2 
PALO  ALTO  CA  94304 

1  OLIN  ORDNANCE 

ATTN  V  MCDONALD  LIBRARY 
P  O  BOX  222 
ST  MARKS  FL  32355-0222 

1  PAUL  GOUGH  ASSOCIATES  INC 

ATTN  P  S  GOUGH 
1048  SOUTH  STREET 
PORTSMOUTH  NH  03801-5423 

1  HUGHES  AIRCRAFT  COMPANY 

ATTN  T  E  WARD 
PO  BOX  11337 
TUCSON  AZ  85734-1337 
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1  SCIENCE  APPLICATIONS  INC 

ATTN  R  B  EDELMAN 
23146  CUMORAH  CREST 
WOODLAND  HILLS  CA  91364 

3  SRI  INTERNATIONAL 

ATTN  G  SMITH 
D  CROSLEY 
D  GOLDEN 

333  RAVENSWOOD  AVENUE 
MENLO  PARK  CA  94025 

1  STEVENS  INSTITUTE  OF  TECH 

DAVIDSON  LABORATORY 
ATTN  R  MCALEVY  ffl 
HOBOKEN  NJ  07030 

1  SVERDRUP  TECHNOLOGY  INC 

LERC  GROUP 

ATTN  R  J  LOCKE  MS  SVR  2 
2001  AEROSPACE  PARKWAY 
BROOK  PARK  OH  44142 

1  SVERDRUP  TECHNOLOGY  INC 

ATTN  J  DEUR 

2001  AEROSPACE  PARKWAY 
BROOK  PARK  OH  44142 

3  THIOKOL  CORPORATION 

ELKTON  DIVISION 
ATTN  R  BIDDLE 
R  WILLER 
TECH  LIB 
P  O  BOX  241 
ELKTON  MD  21921 

3  THIOKOL  CORPORATION 

WASATCH  DIVISION 
ATTN  S  J  BENNETT 
P  O  BOX  524 

BRIGHAM  CITY  UT  84302 

1  UNITED  TECHNOLOGIES  RSCH  CENTER 

ATTN  A  C  ECKBRETH 
EAST  HARTFORD  CT  06108 

1  UNITED  TECHNOLOGIES  CORP 

CHEMICAL  SYSTEMS  DIVISION 
ATTN  R  R  MILLER 
P  O  BOX  49028 
SAN  JOSE  CA  95161-9028 


1  UNIVERSAL  PROPULSION  COMPANY 

ATTN  H  J  MCSPADDEN 
25401  NORTH  CENTRAL  AVENUE 
PHOENIX  AZ  85027-7837 

1  VERITAY  TECHNOLOGY  INC 

ATTN  E  B  FISHER 
4845  MTLLERSPORT  HIGHWAY 
EAST  AMHERST  NY  14051-0305 

1  FREEDMAN  ASSOCIATES 

ATTN  E  FREEDMAN 
2411  DIANA  ROAD 
BALTIMORE  MD  21209-1525 

3  ALLIANT  TECHS  YSTEMS 

ATTN  C  CANDLAND 
L  OSGOOD 
R  BECKER 
600  SECOND  ST  NE 
HOPKINS  MN  55343 

1  DIRECTOR 

US  ARMY  BENET  LABS 
ATTN  SAM  SOPOK 
AMSTRA  AR  CCB  T 
WATERVLIET  NY  12189 
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36  DIR  USARL 

ATTN:  AMSRL-WT-P,  A  HORST 

AMSRL-WT-PC, 

R  A  FIFER 
GF  ADAMS 
W  R  ANDERSON 
R  A  BEYER 
S  W  BUNTE 
C  F  CHABALOWSKI 
K  P  MCNEILL-BOONSTOPPEL 
A  COHEN 
R  CUMPTON 
R  DANIEL 
D  DEVYNCK 
N  F  FELL 
B  E  FORCH 
J  M  HEIMERL 
A  J  KOTLAR 
MRMANAA 
W  F  MCBRATNEY 
K  L  MCNESBY 
S  V  MEDLIN 
MS  MILLER 
A  W  MIZIOLEK 
S  H  MODIANO 
J  B  MORRIS 
JE  NEWBERRY 
S  A  NEWTON 
R  A  PESCE-RODRIGUEZ 
B  M  RICE 
R  C  SAUSA 
M  A  SCHROEDER 
J  A  VANDERHOFF 
MWENSING 
A  WHREN 
J  M  WIDDER 
C  WILLIAMSON 
AMSRL-CI-CA,  R  PATEL 
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USER  EVALUATION  SHEET/CHANGE  OF  ADDRESS 


This  Laboratory  undertakes  a  continuing  effort  to  improve  the  quality  of  the  reports  it  publishes.  Your  comments/answers 
to  the  items/questions  below  will  aid  us  in  our  efforts. 

1.  ARL  Report  Number  ARL-TR-981 _ Date  of  Report  March  1996 _ 

2.  Date  Report  Received _ 

3.  Does  this  report  satisfy  a  need?  (Comment  on  purpose,  related  project,  or  other  area  of  interest  for  which  the  report 

will  be  used.) _ 


4.  Specifically,  how  is  the  report  being  used?  (Information  source,  design  data,  procedure,  source  of  ideas,  etc.) 


5.  Has  the  information  in  this  report  led  to  any  quantitative  savings  as  far  as  man-hours  or  dollars  saved,  operating  costs 
avoided,  or  efficiencies  achieved,  etc?  If  so,  please  elaborate. _ 


6.  General  Comments.  What  do  you  think  should  be  changed  to  improve  future  reports?  (Indicate  changes  to 
organization,  technical  content,  format,  etc.) _ 


Organization 

CURRENT  Name 

ADDRESS  _ 

Street  or  P.O.  Box  No. 

City,  State,  Zip  Code 

7.  If  indicating  a  Change  of  Address  or  Address  Correction,  please  provide  the  Current  or  Correct  address  above  and  the 
Old  or  Incorrect  address  below. 


Organization 


OLD  Name 

ADDRESS  _ 

Street  or  P.O.  Box  No. 


City,  State,  Zip  Code 


(Remove  this  sheet,  fold  as  indicated,  tape  closed,  and  mail.) 
(DO  NOT  STAPLE) 
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POSTAGE  WILL  BE  PAID  BY  ADDRESSEE 


DIRECTOR 

U.S.  ARMY  RESEARCH  LABORATORY 
ATTN:  AMSRL-WT-PC 

ABERDEEN  PROVING  GROUND,  MD  21005-5066 
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